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ABSTRACT 

A thermodynamically admissible theory of viscoplasticity with two internal variables (a 
back stress and a drag strength) is presented. Six material functions characterize a specific 
viscoplastic model. In the pursuit of compromise between accuracy and simplicity, a model 
is developed that is a hybrid of two existing viscoplastic models. A limited number of appli- 
cations of the model to Al, Cu, and Ni are presented. A novel implicit integration method is 
also discussed, and applications are made to obtain solutions using this viscoplastic model. 


NOMENCLATURE 

material constants 

a,A,A\b,c,D 0 ,E,h,H l ,H 2 ,k,ni,n,n , ,Q,T lt a,v 

magnitudes 

PI, IS I, IS-* I, |S-X|,|Xm, W v 

internal variables 
direction tensors 

dij * » Wy * Xy 

material functions 
l » L , r , R , Z , Z„ , 9 

strains and stresses 

^ij > Sy, £ij, £y, £{j, Olj 

1. INTRODUCTION 

The design of structures exposed to high-temperature environments requires mathematical models capable 
of accurately predicting short-term plastic strains, long-term creep strains, and the interactions between them. 
Multiaxial, cyclic, and nonisothermal loading histories are normal service conditions, not exceptional ones. In 
the pursuit to achieve this objective, a hybrid model is developed using walker’s [1] theoretical structure for 
the evolution of internal state and miller’s [1] methodology for characterizing material functions. The main rea- 
son for doing inelastic analysis is not so much to determine the deformation response of the structure per se, but 
rather to assess its useful service life. Results from inelastic analysis are input to life assessing schemes. These 
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schemes typically require such quantities as: stress range, mean stress, inelastic strain range, and the percentages 
of time-dependent creep and time-independent plasticity that make up the inelastic strain range [3]. It is there- 
fore important that the mathematical model used in high-temperature structural analysis be capable of predicting 
these parameters with reasonable accuracy. 

A wide range of materials are used in structures that are designed to operate in high-temperature environ- 
ments. The designer’s choice of material is usually dictated by cost and by such factors as its ability to: carry 
load, conduct heat, resist crack initiation and growth, and resist oxidation and corrosion. A typical metal may 
have some desirable properties, and some undesirable ones as well. With the technology of metal-matrix fabri- 
cation coming of age, the engineer can now begin to tailor materials that provide optimal characteristics for a 
given application. At NASA, composite matrices of aluminum, copper, nickel, and their alloys are being con- 
sidered for a variety of applications. Along with ceramics, tungsten and tungsten alloys are being considered as 
fiber materials. An objective of this paper is to develop a model for characterizing their viscoplastic behavior. 
These material representations will be used in the near future in a mesomechanics model [4,5] to predict the 
viscoplastic response of metal matrix composites. 

The need for reliable and expedient algorithms to numerically integrate viscoplastic models is of 
paramount importance. Such models are composed of systems of stiff, first-order, linear, differential equations 
whose coefficients are coupled and strongly nonlinear, hence, they are inherently difficult to integrate. The 
implicit integration algorithm of walker [6,7], which is presented in appendix A and extended there to include 
the second-order term in its difference formulation, is applied to our viscoplastic model. 

2. VISCOPLASTIC THEORY 

The structure of our viscoplastic theory is a special case of a much more general structure that freed & 
chaboche [8] have shown to be physically motivated and thermodynamically admissible. The theory admits two 
fundamental internal variables; they are: the (positive scalar- valued) drag strength D > 0, and the (deviatoric 
tensor- valued) back stress B^. The drag strength accounts for isotropic hardening effects. The back stress 
accounts for kinematic (or flow-induced anisotropic) hardening effects, and is considered to be the sum of two 
separately evolving constituents. One constituent represents short-range back stress effects, while the other 
constituent B t j represents long-range back stress effects [9,10]. These internal variables are considered to 
evolve phenomenologically through competitive processes associated with strain hardening, strain-induced 
dynamic recovery, and time-induced thermal recovery. 

The strain e, ; is taken to be the sum of an elastic (or thermodynamically reversible) strain £,' and an 
inelastic (or thermodynamically irreversible) strain tfj, i.e., 

Ey = £y + (1) 

with no inelastic strain occurring in the stress-free virgin state. Small strains, displacements, and rotations are 
considered to make up the deformation of a material element. 

The constitutive equations of the theory are given by: 

l-j-V v 

£.7 = —jT a ij ~ £ CT « 5 v + a (T-T 0 )8ij (2a) 

¥=1 3 -yt=i 

D=D 0 + hS (2c) 

where £ is the elastic modulus, v is the Poisson ratio, a is the coefficient of thermal expansion, T-T 0 is 
the change in temperature, is the stress, and 8,y is the Kronecker delta. Here p$ and 5 > 0 are the ther- 
modynamic displacements conjugate to the thermodynamic forces B,J and D-D 0 > 0, with H y > 0 and 
h > 0 denoting their hardening coefficients, and D 0 > 0 defining the annealed value of drag strength. Repeated 
Latin indices are summed over from 1 to 3 in the usual manner. 

The evolution equations of the theory are given by: 

e{j = QZ ntj = |e p | 

P/7 = d., - 0 R Xij 
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(3a) 

(3b) 


(3c) 


P.7 = 

f N 

6= ttV'T ^ ~ Qr (3d) 

U—Vq l 

k . 

which, when combined with the constitutive equations (2b and 2c), give rise to the governing viscoplastic equa- 
tions of our theory; they are: 

ifj = 9 Z «y = |e^| n u (4a) 

|e p | dij+QRxij + 


B i} = |(//,+//2) itj - |Wi -f 1 



(4b) 

(4c) 


where 


X tJ = fly -~H 2 z* = B,j 


V* =(5/y -By)eg = |S-B||e p | 
with directions (unit or projection) defined by: 

_ 3 5y — By _ 3 5y— Xy _ 3 

n ‘ j ~ 2 |S-B| ’ UiJ ~ 2 |S-X| ’ Xij ~ Ypn 

, . ti u X u 

dq — (1 ~ a) Xij + a tty 

and with norms (or magnitudes) defined by: 

|e p | = V(2/3)e5ef ; , [/| = V(3/2>Vy 


(5a) 

(5b) 

(5c) 

(5d) 

(5e) 


where /ye [Sy fly ,Sij-Bij Jtij JSij-Xy } . Sy = <Jy - Ofct8y/3 is the deviatoric stress, and Sy— By is the viscous 
stress that drives inelastic flow. This choice for the directions and norms scales the theory for tension. The 
notation Xy=fly is used for easier visual distinction between Xy and By. The derivatives of the hardening 
coefficients /i, H\ and H 2 represent their change with temperature, if any. 

The material functions of the theory are: 6(T) > 0 is the thermal diffusivity, Z (|S -B \!D ) = |e p |/0 > 0 is 
the Zener [11] parameter, L > 0 and / >0 are the limiting states associated with dynamic recovery, and 
R > 0 and r > 0 are the thermal recovery functions. At this point in the development, the limit and thermal- 
recovery functions are general functions of state. Specific forms for the functions Q,Z,l,L,r, and R will 
vary from model to model. There is one material constant in these evolution equations; it is the nonproportional- 
ity parameter ae[0,l] found in eqn. (5d). This parameter provides added flexibility in modeling nonpropor- 
tional kinematic behavior; it has no effect on proportional behavior. What it does (more or less) is proportion 
the nonlinear strain-induced characteristics of the back stress between the mechanisms of hardening and recovery 
[8]. The value of this parameter has a profound impact on predicted response under nonproportional ratchetting 
conditions [12]. 

Notice that the drag strength D evolves with W V /(D-D 0 ). It cannot evolve with either |e p | or 
(the choices in existing models), aid satisfy the thermodynamic constraint of intrinsic dissipation throughout the 
entire domain of state space [8] defined by the set [T, < 7 y , D , B y). T herefore, in the absence of recovery, the 
drag strength parabolically hardens with viscous work, i.e., D « vW*. Parabolic hardening of the drag strength 
was introduced phenomenologically by MILLER & SHERBY [13] into their model; here it is introduced as a direct 
consequence of the thermodynamics. Because D-D 0 is in the denominator of this forcing function, it appears 
to be singular in the virgin state, D =D 0 . But actually the forcing function is zero valued there, because the 
numerator, W v , is zero valued in the virgin state, too. 
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The basic form of the evolution equations (3b and 3c) for the kinematic variables p,} and was pro- 
posed by walker [1], and is a compromise between accuracy and simplicity. By assuming that the long-range 
back stress B$ does not thermally recover, which is physically incorrect (a compromise), it can be integrated 
exactly and combined with the short-range back stress B,j = X,j. The result is a single governing equation (4b) 
for the back stress (instead of two) which has both short- and long-range features providing improved accuracy. 
The direction for dynamic recovery d tJ of the back stress, as established in eqn. (5d), presents a refinement to 
the original theory. Dynamic recovery is taken to be proportioned by a factor a between directions opposing 
the short-range back stress and the stress difference Sy-Xy. The scaling factor «yXy/|X| in eqn. (5d) is 
introduced to make the model insensitive to the direction of dynamic recovery under proportional loading condi- 
tions (i.e., the predicted response is the same for all values of the parameter ae [0,1]). The model only becomes 
sensitive to the direction of dynamic recovery (i.e., the value of a) under nonproportional loading conditions. 
This scaling factor is also required for the theory to be thermodynamically admissible [8]. 


3. VISCOPLASTIC MODEL 


Our model is predicated upon the assumption that the magnitude of inelastic strain rate (in particular, Z ss ) 
can be represented as a sum of two power functions in stress (cf. appendix B). One to govern at the higher- 
temperatures/lower-stresses, while the other to govern at the lower-temperatures/higher-stresses. An analogous 
expression for Z has recently been proposed by chaboche [14] to model the transient strain-rate sensitivity of 
a stainless steel. This parameter is a generalization of walker’s [1] power-law model, and a simplification 
(especially from a numerical integration viewpoint) of miller’s [2] hyperbolic sine model, i.e., a compromise 
between accuracy and simplicity. We motivate its choice with the steady-state creep data of fig. 1, where 
Z« = |e c |/0. By applying Miller’s hypothesis (cf. appendix B) to our relationship for the stress dependence of 
inelastic flow at steady state, we are able to obtain explicit relationships for the transient Zener parameter Z, 
and both of the thermal recovery functions r and R . Details on the derivation of the material functions, and a 
procedure for determining values for their associated constants, are presented in appendix B. 


The material functions of our viscoplastic model are: 
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(6a) 


(6b) 

(6c) 

(6d) 

(6e) 

(6f) 


where eqn. (6e) implies no dynamic recovery of the drag strength. (A particular limit function l is under con- 
sideration to account for monotonic/cyclic interactions, which are not incorporated in the current model, i.e., 
I = oo.) The material constants for this model are comprised of: the elastic constants: E , a, and v; the thermal 
activation energy and transition temperature: Q and T , ; the steady-state constants: A, A', n, and n the 
constant-structure constant: m\ the threshold constants: b and c; the transient constants: h,H x , and H 2 ‘, the 
nonproportional constant: a; and the initial values: D 0 and D(t = 0). The parameter k is the Boltzmann con- 
stant, and H(D-D 0 ) is the Heaviside function. Values for these constants are given in table 1 for aluminum, 
copper, nickel, and tungsten. If one chooses, a reference temperature could be introduced (thereby setting 
0(7’„/)= 1) which would, in effect, scale the drag strength. 
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This model has several deficiencies that the authors are aware of, and are working on. It does not account 
for monotonic/cyclic interaction effects. It overpredicts normal ratchetting behavior. And it does not account for 
inverse strain rate sensitivity effects exhibited by many alloys. These deficiencies are common amongst most 
viscoplastic models found in the literature. 


4. NUMERIC INTEGRATION 

One can combine the constitutive equations (2) and the evolution equations (3) of our viscoplastic theory 
and obtain the following system of stiff, linear, first-order, difference equations: 

AL,j + E, y Aa = Axij (7a) 

A Xij + Xij A b = Ay iy (7b) 

AK + K Ac = Az (7c) 


with 


Sij = £>y + Btj 


Bij = Xij + y #2 e y 


e p - 

c ‘j 


e ij 


2p 

D = Do + '&K 

where the coefficients are nonlinear and coupled, and are defined by: 

A a = 


= 3 ^ |Ae p | - — Au 
\S-B\ 1 p * 


Ab =Hi 


|Ae p | + 0 — At. 

L 1 1 \X\ | //, 


AH i 


Ac =2 


IJ^BI |Ae p , + 0-i- 

/ D-D o D-D 0 


At 




and where the nonhomogeneous terms are also nonlinear and coupled, and are defined by: 


B 


A = 2p A - ABij + — ^ L A|i 


4 


Ay./ = f 


- a 


UklXkl 


|Ae p | 


Az = ft |S-B | |Ae p | 


(8a) 

(8b) 

(8c) 

(8d) 

(9a) 

(9b) 

(9c) 

00a) 

(10b) 

(10c) 


Here A|i, Ah , and AH t represent the incremental change in these moduli due to an incremental change AT in 
temperature. This system of difference equations characterizes the deviatoric response of the material, while the 
constitutive equation a u = 3k[£ u - a(T-T 0 )8 u ] characterizes the hydrostatic response of the material. Here 
(j. = £/2(l+v) is the elastic shear modulus, and k = £/3(l-2v) is the elastic bulk modulus. 

The system of difference equations (7) can be solved by classical methods, such as the explicit methods of 
Euler and Runge-Kutta. Such methods require numerous time steps in order to maintain numeric stability and 
accuracy because the coefficients and nonhomogeneous terms make this system of equations mathematically 
stiff. This is not a detriment when numerous data pairs are needed to plot smooth accurate response curves. 
However a new method, the uniformly-valid implicit integration scheme of walker [6], may be preferred for 
finite element analyses where computational times are large, and the ability to create smooth response curves at 
the nodal points is not the typical required output [7]. An advantage of this method is that it is not constrained 
by time step size in order to maintain numeric stability and accuracy since the scheme is uniformly valid. 
Walker’s method transforms the system of difference equations (7) into recursive integral equations, which are 
then solved using an asymptotic expansion technique (cf. appendix A). The resulting integration scheme for the 


5 


general viscoplastic theory presented in eqn. (4) is given by the following system of stiff, coupled, implicit, 
recursive relations: 


XijO+At) = Sy(0 + 


Axij(t+At) 

Aa(t+At) 


Axj/J+At) - AXjjjt) 
Aa (t +At) 

X ij (t+At) = X ij (t)e 

Ay jj(t+ At) - Ay tJ (t ) 
Ab(t+At) 


V -Aa(t+At) _ 


1 


\ - e - Aa(,+A,) j + 

1 _ g-Aa(l+A<)] 


-Afr((+AI) , Ay ij (t+At) 
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1 _ g-Ai(t+A<) 


e 


-Ab(t+At) _ 


l 


Ab(t+At) 


1 _ g-A*(l+A() 


K(t+At) = K(t) \ 1 - e _A£(,+A,) ] + 

Ac(t+At) ^ J 


Az(t+At) - Az(t ) 
Ac ( t+At ) 


-Ac (t+At) 


l 

Ac (t+At) 


l _ g-Ac(<+A<) 


(Ha) 


(lib) 


(He) 


where t > 0 and where At > 0 may be arbitrarily large or small. These equations must be solved simultane- 
ously using a technique such as Newton-Raphson iteration, where Aa, Ab, and Ac are the iterated variables 
[6,7]. For viscoplastic models that have a back stress and a drag strength (or yield strength) which evolve, 
Walker’s implicit integration scheme has only three variables to iterate over; whereas, the classical implicit 
integration schemes (e.g. backward Euler difference) have thirteen variables to iterate over (i.e., six for inelastic 
strain, six for back stress, and one for drag strength). 


5. APPLICATIONS 

A few applications for aluminum, copper, and nickel (data are insufficient to include tungsten) are 
presented to demonstrate some of the capabilities and limitations of the model, and to show feasibility of the 
Walker integration method. 

The model’s capability in predicting the monotonic behavior of aluminum as a function of temperature is 
demonstrated in fig. 2. The model produces acceptable predictions above approximately .2 T h (or about - 
100°C, where T h is the material’s homologous temperature). This is an improvement over a single power-law 
for the Zener parameter with an exponent of n~ 5, which produces acceptable predictions above approximately 
.5 T h (or about 200°C for aluminum). In applications requiring modeling at both elevated temperatures and 
temperatures below .2 T h , the hyperbolic Zener parameter of miller [2] is recommended. Numeric integration 
using the technique of Walker described in the previous section and the second-order method of Runge-Kutta 
with automatic time-step control was done here with good agreement between the two methods (less than 1 MPa 
difference, with a maximum error of 20% which occurred at the highest temperatures where the stress is on the 
order of 5 MPa in fig. 2). The method of Walker was found to be superior to that of Runge-Kutta in both accu- 
racy and computation time. This was especially true when the model became mathematically stiff (i.e., in the 
regions where hardening and recovery mechanisms are of comparable magnitude). Here the difference in cpu 
times exceeded an order in magnitude, and the Runge-Kutta solutions oscillated. 

A thermomechanical prediction of the model (in this case, for nickel) is presented in fig. 3a. The model 
produces the correct shape of the hysteresis loop with two exceptions. Of course, it cannot model the serrated 
yielding (or strain avalanching, [28]) observed in the experimental response. And second, the model presently 
has no provision to account for monotonic/cyclic interactions, and therefore it continues to isotropically harden 
beyond that which is experimentally observed. This deficiency is being worked on. 

Figure 3b presents a comparison between predictions obtained using the methods of Walker (with several 
different time steps) and second-order Runge-Kutta (with automatic time stepping). The Walker analysis with 
the largest time step size (50°C / step) ran about 2.5 times faster than the Runge-Kutta analysis, the one with the 
intermediate time step (approximately 15°C / step) ran about 1.5 times faster, while the one with the smallest 
time step (1.5°C / step) ran about 3 times slower. The greatest variation in the three Walker solutions occurs at 
the knees of the hysteresis curve. A benefit of the Walker (implicit integral) method is that this error is not 
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propagated, as all three solutions agree at the turnaround points. This is not true of explicit differential methods, 
such as Runge-Kutta, where an error initiated is also propagated. The effect of a slowly accumulated propagat- 
ing error is found in the Runge-Kutta solution, and can be observed by looking at successive peaks in fig. 3b. 

The model’s capability in predicting nonproportional material behavior (for copper) is demonstrated in fig. 
4. Here the influence of the nonproportionality parameter a on predicted response is shown. The limiting 
case of a = 0 represents the nonlinear kinematic hardening rule of Armstrong and Frederick [9,14], while the 
other extreme of a = 1 denotes a nonlinear Prager-like hardening rule [8,12]. The prediction for a = 1 pro- 
duces an unwanted ratchetting down the axial stress axis; a result also observed by lamba & sidebottom [29] 
in plasticity analyses using the Prager and Ziegler hardening rules. There is little difference in the predicted 
response over a wide range in the parameter a (a = 0 to a ~ 0.9 in this experiment). Therefore a more 
severe nonproportional loading history (such as advocated in ref. 12) is necessary to characterize the value of the 
material constant a . This figure illustrates another limitation of the model. Our model, like most, is formu- 
lated with a von Mises type norm, whereas copper behaves like a Tresca material [29]. That is why the 
predicted axial stress response is too soft, as the value of the drag strength was correlated with the shear 
response in this figure. 

6. CONCLUSIONS 

Within the framework of a thermodynamically admissible theory of viscoplasticity with an evolving back 
stress and drag strength, a model is presented that is a hybrid of the walker [1] and miller [2] models. The 
model is for pure metals, and has been applied to aluminum, copper, and nickel. Tungsten has been partially 
characterized, using rules of thumb instead of experimental data (they do not exist) to establish values for 
roughly half of the constants. Both strengths and weaknesses of the model are discussed. The implicit integra- 
tion scheme of walker [6,7] has been further developed and applied to our model, and was found to be superior 
to the second-order method of Runge-Kutta. A detailed numerical study of the Walker integration technique is 
currently underway. 
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APPENDIX A 


WALKER IMPLICIT INTEGRATION METHOD 


A recursive integration relation for a linear, first-order, differential equation is given. The uniformly-valid 
asymptotic expansion of walker [6] is applied to this integral equation. In the past, this method has been 
applied to specific viscoplastic models retaining only the first term in the expansion [6,7]. Here we apply it to a 
general differential equation so that its structure is more fully revealed, and it is extended to include the second 
term in the expansion. 

Consider the linear, first-order, differential equation 


X +u X = v 


(Al) 


where u and v are the parameters, which may be nonlinear functions of X, and X is the variable to be 
solved for. One can integrate this differential equation over the time interval [0,f+Af] and obtain 


t+At 


X(t+ Af) = X>~“ ( ' +A °+ J exp 


$=0 

t+At 


t+At 


-J 






d % 


(A2) 


= ^ 0 g-“( ,+A ») + f e ~lu(l+Al) - u(5)] dv(%) 

where X 0 = X (0). The integral in this equation is the nonhomogeneous contribution to the solution of the 
differential equation (Al). This integral equation can be cast into a recursive relation by dividing the interval of 
integration into two parts, viz., 

X(t+At) = X 0 e- u(,+A,) + (*-[«('+*> -«©1 + 




t+At 


+ f e -lu(t+ At) - u(£)l dv(k) 

£=t <% ' ' 


(A3) 


Substituting the identity 1 = e into the first integral of this equation results in 

X(t+At) = AV~“ (,+a ° + e ~ lu(t+A,) - “<'» f e -[«(0 + 


t+At 


r e -[u(,+At) - *(6)1 dv(5) ^ 

5=f 

which simplifies to the desired recursive integral equation 

X(t+At) = X(t)e~ lu(,+A,) ~ u(,)] + f g -["(M-A<)-«( 5 )] j y(0 


k=< 


(A4) 


(A5) 


where t > 0 and u (0) = X (0) = 0. The derivation of this relation made use of eqn. (A2) when integrated over 
the time interval [0,r]. 

The recursive relation given by eqn. (A5) for the differential relation in eqn. (Al) becomes practical only 
when a solution strategy exists for evaluating the integral 

/(Ar)= J (A6) 

When u(t) is a monotonically increasing function of /, this integral possesses an evanescent memory of the 
forcing function 3v(^)/3|, causing most of the contribution to the integral to come from the region near the 
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upper limit ^ ~ f+Ar, where the exponential term is approximately equal to its maximum value of 1. Using the 
method of Laplace, walker [6] derived a uniformly-valid asymptotic expansion for this integral equation by 
first expanding the arguments of the exponential and derivative into a Taylor series, and then integrating. In par- 
ticular, by writing u© = a [f+Af - (r+Ar-£)] and expanding it by Taylor’s theorem to give 

u© = u(t+At) - (r+A/©) + • ■ • , the integral equation (A6) becomes 


t+At 


I (At) = f e -[«+A»-§«- •••]. 


dvfe) 


d% 


which can be rewritten in equivalent form as 


I (At) = 2) 


dz 


(A7) 


(A8) 


z-=0 ^ 

by introducing the change in variable z = t+At-Z,. If one now expands the argument of the derivative in this 
integrand using Taylor’s theorem to give v(t+At-z) = v(t+At) - z + • • • , then this derivative can 

be represented by the series **v(t+At z) _ _ 1 ) = - y + zv - • • • , and the integral in eqn. (A8) 

can be expressed as 




I(At)= Je~<* >(v - zv + 


)dz 


(A9) 


z=0 


Integrating this relationship gives a uniformly valid asymptotic expansion for the integral in eqn. (A6); v/z.. 


/(A/) = -fl - + ~ At e-“ A ‘ - Lfi 

u > ' i) i) ^ J 


(A10) 


where the third and higher-order terms in the expansion have been neglected. Here we recall that the deriva- 
tives of u and v are evaluated at time /+ At. 

A uniformly-valid, implicit, numeric, integration method providing a solution strategy to the recursive 
integral equation (A5) for the linear, first-order, differential equation (Al) is therefore given by the difference 
relation 

X(t+At) =X(t)e- Au(,+A,) + [l " «" A “ ( ' +A °) + 


Av(t+At) - Av (t) 


Au(t+At) 


-Au(t+At) 1 f i _ --Au(l+A«)| 

Au(t+At)( J 


(All) 


where Ax(t+At) = x(t+At) - x(t) ~ x(t+At)At , Ax(t+At) =x(t+At) - x(l) = x.(t+At)At, and hence 
Ax(t+At) - Ax(t) ~ [x(t+At) - x(t)]At ~ x(t+At)(At) 2 for a constant time step size At. This implicit 
integration scheme must be solved iteratively, since Au(t+At) and Av (t +At ) may be nonlinear functions of 
X(t+At). 

The implicit integration method of eqn. (All) is valid for an arbitrary time step size At > 0, be it large 
or small, even when only the first term in the expansion is retained, i.e., it is uniformly valid. This follows 
because whenever A«»l, or equivalently An X»AX, eqn. (All) reduces to 


f ('“) = 77TT 

Au(t+At) 

and whenever And, or equivalently A u X cAX , it reduces to 

X(t+At) ~X(t) + Av(t+At)~X(t) + AX(t+At) 
when keeping only the first term in the asymptotic expansion, or to 

X(r+Ar) ~X(t) + Vz[Av(t+At) + Av(f)] 

= X(t) + V2[AX(t+At) + AX (01 


(A12) 


(A 13) 


(A 14) 


10 


when keeping the first two terms in the asymptotic expansion, as in eqn. (All). Equation (A12) is the asymp- 
totic solution to eqn. (Al) whenever Am X>AX, and eqns. (A13 and A14) are valid whenever Am XcAX. In 
fact, eqn. (A13) is the Euler backward difference method, while eqn. (A14) is the trapezoidal difference method. 
Because keeping the second term of the asymptotic expansion comes with negligible computational expense per 
iteration (one must only retain Av(0 from the previous time step), and because numerical experiments show 
that this reduces the number of iterative steps required to attain convergence, it is therefore recommended that 
the second term in the asymptotic expansion be retained in numeric algorithms. 

Viscoplastic models are composed of systems of stiff, linear, first-order, differential equations with coupled 
nonlinear parameters. An algorithm for the numeric solution of these models using the implicit integration 
scheme of eqn. (All) (where only the first term in the asymptotic expansion was kept) and the iterative tech- 
nique of Newton-Raphson is given in ref. 6. An earlier algorithm of walker’s [1] (where the derivatives of u 
and v are evaluated at time t instead of time t+At) is not as efficient. 

Even though this integration scheme is stable and accurate for time steps of arbitrary size from the numer- 
ics viewpoint, it will still be necessary to moderate the time step size from a physics viewpoint in order to prop- 
erly trace the loading history, especially when it is nonproportional such that the principal axes of tensor-valued 
variables rotate in a material reference frame. The evolutionary equations in unified viscoplastic formulations 
introduce an evanescent strain memory effect (evanescent along the inelastic strain path). This effect is con- 
sistent with il’YUSHIN’s [30] delay-trace hypothesis, that a material’s response depends not on the whole of the 
previous strain path, but only on the most recent part of it. This physical fact is assimilated in the integration 
algorithm by expanding the exponential integrand in eqn. (A6) about the current state. 

APPENDIX B 

MODEL DEVELOPMENT AND CHARACTERIZATION 

Here we derive the material functions given in eqn. (6) that are unique to this model, and present a metho- 
dology for determining values for the material constants of the model. 


MILLER’S HYPOTHESIS 


The Zener parameter Z in the flow equation (3a) is a temperature normalized function describing the 
magnitude of inelastic strain rate, i.e ., Z = |e p [/0. At steady state, one must be able to express Z as only a 
function of stress, since the internal variables do not evolve at steady state. Therefore, one can propose two dis- 
tinct Zener parameters, each differing in argument and function; one for steady states, and the other for transient 
states. If one defines these arguments as 




’\S\ 

7 — 7 

\S -B\' 

A 

K d 

, — Id 

D 


then the hypothesis of Miller [2] given by 


Z = Z w 




f 15 -Bl' 


D 

- j 



(Bl) 

(B2) 


establishes the functional dependence of the transient Zener parameter Z (the harder of the two to characterize) 
in terms of the steady-state Zener parameter Z w . Here A > 0 is a material constant introduced to normalize 
the units, and m > 1 is the constant-structure exponent. 

The advantages of this hypothesis go beyond just providing the functional form for the transient Zener 
parameter. One can now equate arguments between the Z of eqn. (B2) and the Z„ of eqn. (Bl) under steady 
state conditions. Assuming that the steady-state value for back stress is proportional to the value of applied 
stress, i.e.. 


then it immediately follows that 


\B\ss=b |S | 


=A(l-b) 



m-1 


m 


(B3) 

(B4) 


11 


independent of the choice of functional form for Z„. Actually, Miller proposed eqn. (B4) and then derived eqn. 
(B3), where the proposition of eqn. (B4) was based on data from warm-worked material. These relationships 
play an important role in characterizing the thermal recovery functions. Note that m = 1 implies 
= A (1-7?) = constant, but in the model of eqn. (6) m > 1 because of the recovery function r. 


MATERIAL FUNCTIONS 


A steady-state Zener parameter that is the sum of two power-law functions, i.e.. 


\s\ 

n 

A 

n' 

\S\ 

A 

t J 


A’ 

i. J 


A 

t. j 


(B5) 


is proposed, and we motivate this choice with the steady-state creep data of fig. 1. Applying Miller’s 
hypothesis, eqn. (B2), to this form for Z ra results in the transient Zener parameter Z given in eqn. (6b). The 
function for thermal diffusivity, eqn. (6a), was derived by MILLER [2], 

The limit function L for the short-range back stress controls the degree of curvature exhibited at the 
knee of the stress-strain curve. Equation (6c) assumes that L is proportional to the value of the drag strength. 
This relationship is in accordance with the physical observation that a material’s capacity to sustain an internal 
stress (of which L is the measure) is directly related to the density and distribution of its dislocation substruc- 
tures (of which D is the measure). 

By definition B t j =0 at steady state, and therefore from the constitutive equation (2b) one obtains 

= =0 (B6) 


from which one derives the thermal recovery function R for the back stress. Here it is considered that T = 0 
at steady state, which implies that H X = H 2 = 0 at steady state. Combining this relationship with the evolution 
equations (3) for py and P -j, and using the fact that X j; , By, and Sy are coaxial with each other at steady 
state (a direct consequence of the chosen evolution equations), leads to 


R = 


H } +H 2 

h7 



(B7) 


which when combined with eqns. (B3 and B5) results in the expression for the thermal recovery function R 
given in eqn. (6d). 

At steady state, D = 0 by definition, and therefore the evolution equation (3d) for drag strength reduces 
to 

"Z» (B8) 


0 r = 


D \S-Bj 


D-D c 


le'l 


= 0 


D 


D-Dr 


IS I 


where eqns. (B3 and B4) are used to reduce the coefficient |S -B |/D in the middle expression. Combining 
eqns. (B4, B5, and B8) results in the thermal recovery function r for drag strength presented in eqn. (60- 


MATERIAL CONSTANTS 

For a complete material description using the model presented herein, there are: three elastic constants: 
E , a, and v; two thermal constants: Q and T, ; four steady state constants: A , A ', n , and n the constant- 
structure constant: m; two threshold constants: b and c; three transient constants: h,H u and H 2 , the nonpro- 
portionality constant: a\ and the initial value: D 0 . Their values are given in table 1 for aluminum, copper, 
nickel, and tungsten. There are seventeen independent material constants in all, which to a good approximation 
are temperature insensitive. There are a number of rules of thumb (or estimates) given below that one can use to 
assign values to the various constants if data are unavailable. 

The activation energy for self diffusion Q > 0 and the melting temperature T m > 0 are typical hand- 
book data. A useful estimate is that the transition temperature T, > 0 associated with the activation energy for 
inelastic flow is given by T, ~ 0.6 T m with actual values ranging between 0.5T m and 0.87 m [31], 
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With the thermal diffusivity now characterized, one can determine the steady-state Zener parameter 
Z„ = le^U/Q from creep data and plot it against its associated flow stress, as done in fig. 1. The lines in this 
figure represent the fit of eqn. (B5) to the data. Values for the constants A > 0, A ' > 0, n >0, and n ' > 0 are 
readily obtained by eye, or by the method of least squares. Since creep data are probably the most abundant of 
all high-temperature experimental data, one should have no difficulty in characterizing the Norton creep con- 
stants A and n (n ~ 5, typically). High-stress creep data, however, are sometimes scarce or unavailable. If 
this is the case, one can use the estimates given by: n' ~ 3.5 n and Z co = (A'/A) nn ' ,( - n '~ n) = 10*° j - 1 , with 
values ranging between 3 n to 4n for n ' and 10 9 to 10 12 s _1 for the Zener parameter Z co at cross over (cf. 
fig. 1). 

The remaining steady state constant, the Bauschinger coefficient be (0,1) defined in eqn. (B3), can be 
obtained from a tensile test where a specimen is loaded to its saturation (or steady state) strength o sat , and then 
unloaded sufficient enough to induce some reversed plasticity. One then obtains a value for the Bauschinger 
coefficient from the relationship b ~ (a saJ + )l2a sal , where is the reversed yield strength. A word of 
caution. A small offset strain should be used to establish a n to ensure a realistic value for b . Values for b 
may range from 0.2 to 0.7, with a typical value being 0.4. The larger the value for b , the greater the predicted 
Bauschinger effect. 

There is one more material constant that one can ascertain from creep testing. It is the constant structure 
exponent m > 1 which can be determined from stress-drop tests. Fortunately, this parameter is nearly insensi- 
tive to material selection. Typical values for the constant structure exponent are: m ~ 1.8 for pure metals, and 
m ~ 1.5 for alloys [13]. 

The next step is to obtain those constants one can from stable stress-strain hysteresis loops, where the iso- 
tropic variable has saturated. Ratios of plastic strain range to elastic strain range of 1:1 to 3:1 generally work 
well. The constant H 2 > 0 represents the tangent modulus in the region of linear strain hardening after going 
around the knee of the curve, with H 2 ~H i/100 serving as a rule of thumb. The constant Hi> 0 defines the 
tangent modulus at the onset of inelastic deformation. Our experience from optimization calculations is that 
Hi~EI5 provides a good estimate for its value. The larger the value for the material constant c, the more 
rounded the knee in the hysteresis loop will be. The method of nonlinear least squares [32] is a useful tool, but 
not a necessary one, for determining c and II \ (and D sm since the drag strength is an unknown constant at this 
time) since there are three constants to be determined simultaneously at this step of material characterization. 

The constants D 0 > 0 and h > 0 are characterized next. The constant h establishes the rate of isotro- 
pic hardening in monotonic tension to large strains. The initial value D 0 for drag strength is obtained by com- 
paring predicted and experimental values for the yield strength. Some iteration between these two constants 
may be necessary to obtain a proper fit. Again, the method of nonlinear least squares may be helpful, but it is 
certainly not necessary. 

The final material constant to quantify is the nonproportionality parameter a e [0,1], This requires 
nonproportional experimental data. If they do not exist, which is likely to be the case, setting a = 0 implies that 
the theory resembles a two-surface plasticity model with Mroz hardening [33], Further experimental study is 
required before typical values for a can be given. 
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TABLE 1 


Elastic Constants! 

Constant 


Material 

A1 

Cu 

NI 

w 

E o 

GPa 

71 

127 

210 

395 

Ex 

MPa/°C 

-36 

41 

-2.5 

-2.3 

E 2 

MPa/°C 2 

-.011 

-.027 

-.063 

-.027 

V 

- 

.34 

.34 

.31 

.28 

Oo 

1/°C 

25x10"* 

16x10"* 

13x10"* 

4.4x10"* 


1/°C 2 

_ 

5X10" 9 

- 

.16X10" 9 


E = Eq + E \T + E 2 T a 




$ 

a = Oo + a 1 7’ 

T is in degrees centigrade. 



Inelastic Constants* 

Constant 

Units 

Material 

A1 

Cu 

Ni 

W 

A 

MPa 

.076 

1.2 


.29 

A' 

MPa 

4.4 

18 


26 

b 

- 

.6 

.35 


•4t 

c 

- 

23 

7 

15 


Do 

MPa 

.1 

1.2 

.5 

•3t 

h 

MPa 

.05 

3 

2 

5f 


GPa 

10 

20 

70 


H 2 

MPa 

100 

200 

700 


mf 

- 

1.8 

1.8 

1.8 

1.8 

n 

- 

4 

5.5 

5 

4.5 

n' 

- 

14 

19 

15 

19 

Q 

kJ/mole 

140 

200 

290 

590 

T, 

°K 

560 

820 

1000t 

2600 


Data are not available to characterize the material 

* 

constant a 

, nor is there adequate experience to 


assign values via a rule of thumb. 



t 

Determined by engineering judgment because of a lack 
of data. 
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ZENER PARAMETER/ 1/s 



STRESS, MPa 
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MODEL PREDICTION 

EXPERIMENTAL DATA 
(REF. 27) 



(3) FIRST TWO CYCLES AND SATURATED CYCLE. 


WALKER METHOD FOR 
VARIOUS TIME STEP SIZES, 
DEG C/STEP 

1.5 

□ 15 

• 50 


RUN6E-KUTTA METHOD 



MECHANICAL STRAIN 


.02 


s 


CD 

LU 


.01 


-.01 


1,6 



- 0.2 » 1 » 1 » 

-.008 . 004 0 .004 .008 

AXIAL STRAIN, 

(a) STRAIN CONTROL PROGRAM. PATH SEGMENTS SE- 
QUENCE, 0 - 1 - 2 - ... - 8. EACH PATH 
SEGMENT IS APPLIED IN 60s. 


NONPROPORTIONAL 
MATERIAL CONSTANT, 
a 

0 

.75 

.9 

EXPERIMENT (REF. 29) 



(b) THEORETICAL PREDICTIONS USING INTEGRATION 
METHODS OF WALKER AND RUNGE-KUTTA. 

FIGURE 3. - INPHASE THERMOMECHANICAL RESPONSE 
OF NICKEL. CONTROL IS VIA LINEAR RAMPS IN 
TEMPERATURE AND MECHANICAL STRAIN WITH A 
600s CYCLE TIME. TEMPERATURE RANGE, 400 TO 
700 °C. MECHANICAL STRAIN RANGE, ±0.0026. 


(b) EXPERIMENTAL AND PREDICTED STRESS RESPONSE, 
D = 6.5 MPa IN THESE CALCULATIONS. 

FIGURE 4. - CYCLIC, NONPROPORTIONAL, AXIAL- 
TORSIONAL TEST OF COPPER AT ROOM TEMPERA- 
TURE. THE TEST SPECIMEN WAS CYCLED TO 
ISOTROPIC SATURATION PRIOR TO TESTING. 
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